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Abstract 



By the example of sp^-bonded semiconductors, we illustrate what 3rd-generation muffin- 
tin orbitals (MTOs) are. We demonstrate that they can be downfolded to smaller and 
' smaller basis sets: sp^d}^, sp'^, and bond orbitals. For isolated bands, it is possible to 

generate Wannier functions a priori. Also for bands, which overlap other bands, Wannier- 
like MTOs can be generated a priori. Hence, MTOs have a unique capability for providing 
. chemical understanding. 
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1 Introduction 



G Muffin-tin orbitals (MTOs) have been used for a long time in ab initio, e.g. density-functional 
(DF), calculations of the electronic structure of condensed matter. Over the years, several 
^ > MTO-based methods have been developed. The ultimate aim is to find a generally applicable 
Q ', electronic-structure method which is intelligible, fast, and accurate. 

In order to be intelligible, an electronic-structure method must employ a minimal and flexible 
basis of short-ranged orbitals. As an example, the method should be able to describe the valence 
band and the lower part of the conduction band in sp^-bonded materials using merely four short- 
ed [ ranged s- and p-orbitals per atom and, for insulating phases, using merely occupied orbitals such 
as bond orbitals. Another example is materials with strong electronic correlations. For such 
materials, one must first construct a small, but realistic Hilbert space of many-electron wave 
functions, and this requires an accurate and flexible single-particle basis of atom-centered short- 
ranged orbitals. A small basis of short-ranged orbitals is a prerequisite for a method to be 
intelligible and fast, but it may be a hindrance for its accuracy, because the orbitals of a smaller 
basis tend to be more complicated than those of a larger basis. 

Most other density-functional methods, such as plane-wave pseudopotential, LAPW, PAW, 
and LCAO methods, aim at simulation, and are therefore primarily accurate and robust. But 
they are neither fast nor intelligible in the above-mentioned sense, because they employ basis 
sets with of order hundred functions per atom. With such methods, understanding can therefore 
only be attempted after the calculation, by means of projections onto e.g. Wannier functions in 
case of insulators, charge densities, electron-localization functions (ELFs), partial waves, a.s.o.. 

The so-called 3rd-generation MTO method (Andersen et al 1994, Andersen et al 1998, An- 
dersen et al 2000, Tank and Arcangeli 2000, Andersen and Saha-Dasgupta 2000) should come 
close to what we have been aiming for. In the present paper we shall explain what 3rd-generation 
MTOs are and what they achieve. Emphasis will be on the so-called downfolding and energy- 
mesh features which enable MTO bases to be small, flexible, and accurate, as we shall demon- 
strate by exposing them to the above-mentioned sp^-test. From the result, the idea emerges, 
that for band insulators, an MTO basis can be designed a priori to span the Hilbert space of the 



occupied states only. That is, there is one, and only one, such MTO per electron. To get this 
count right, one may associate each orbital with a nominal electron (or pair), and leave it to the 
method to shape the orbitals in such a way that the basis set becomes complete for the occupied 
states. This can be done because MTOs are selective in energy, in the sense that the MTOs 
of order N (NMTOs) are shaped in such a way that the NMTO basis set solves Schrodinger's 
equation exactly for N+1 single-particle energies, which in the present case must be chosen in 
such a way that they span the valence band. This ability to generate Wannier functions directly 
in real space, should be useful for ab initio molecular-dynamics simulations. NMTOs may also 
prove useful for designing many-electron wave functions, which describe correlated electron sys- 
tems in a realistic way. The description of spin and orbital ordering is a trivial example. Also 
for the conduction bands of metals, Wannier-like, low-energy MTOs can be designed a priori. 
This has been demonstrated in several cases (Miiller et al 1998, Sarma and Saha-Dasgupta 2000, 
Korotin 2000, Valenti et al 2001, Dasgupta et al 2002), most recently for the hole-doped cuprate 
high-temperature superconductors, where the material-dependent trend of the hopping integrals 
and their correlation with the maximum Tc was discovered (Pavaxini et al 2001, Dasgupta et at). 

For a description of how we expand the charge density locally in such a way that Poisson's 
equation can be solved and the total energy and forces can be evaluated fast and accurately, 
we refer to previous (Andersen et al 2000, Tank and Arcangeli 2000) and coming publications 
(Arcangeli and Andersen, Savrasov and Andersen). This part of the 3rd-generation computer 
code is still under construction. 

2 Screened spherical waves, kinked partial waves, and muffin- 
tin orbitals 

The 3rd-generation MTO formalism is the multiple-scattering - or KKR (Korringa 1947, Kohn 
and Rostoker 1954)- formalism for finding the solutions, (r) , of Schrodinger's equation for 
an electron in a muffin-tin potential, V (r) = vr (rji) , with the following three extensions: 

1. : The KKR formalism is proved to hold, not only for superpositions of spherically-symmetric, 
non-overlapping potential wells, vr (r^) , but also to leading order in the potential- ouerZap (An- 
dersen et al 1992). Here, and in the following, rR= |r — R| , and R are the sites which wc label 
by R. The potential, vr (r) , is taken to vanish outside a radius, sr, which should not exceed 
1.6 times the radius of touching spheres, i.e.: sr + sr' < 1.6 |R — R'| for any pair of sites, R 
and R'. 

2. : Exact screening transformations of the spherical waves, n/ [nrR) Yl (r^) , are introduced 
in order to reduce the spatial range and the energy dependence (k^ = s) of the wave-equation 
solutions, ipRL (e, r) (Andersen and Jepsen 1984, Andersen et al 1992, Zeller et al 1995). Here, 
and in the following, L = lm labels the spherical- (or cubic-) harmonic's character. 

3. : Energy- independent MTO basis sets are derived which span the solutions (r) with en- 
ergies £i of Schrodinger's equation to within errors proportional to (sj — eo) [si — ei) .. (si — e^r) , 
where eo, ei, cat is a chosen energy mesh with N+l points (Andersen et al 2000, Andersen and 
Saha-Dasgupta 2000) . Such an energy- independent set of Nth-order MTOs is called an NMTO 
set. By virtue of the variational principle, the errors of the energies £j will be proportional to 
{£i - eof {ei - eif .. (e^ - ejv)^ . 

At the top of figure 1 we show the LDA energy bands (k) of Si in the diamond structure, 
calculated with the basis set of Si-centered s-, p-, and d-MTOs, i.e. with 9 orbitals/atom, for 
the 3-point energy mesh eg, ei, 62 indicated on the right-hand side. These bands have meV- 
accuracy for the MT-potential, which in the present case was the standard all-electron DF-LDA 
atomic-spheres potential. Since three energy points were used, the MTOs are of order N=2, 
that is, they are quadratic MTOs, so-called QMTOs. 



The QMTO, Xpin (r) , pointing along [111] from one Si to its nearest neighbor, is shown 
in the (2ll)-plane by the first contour plot. This orbital is localized and smooth, with a few 
"orthogonality wiggles" at the nearest neighbor. The remaining three contour plots show major 
constituents of this pm-QMTO: The pm-kinked partial wave (KPW), (ppiu (e, r) , at the central 




Figure 1: Band structure of Si calculated with the Si spd-QMTO basis set corresponding 
to the energy mesh shown on the right-hand side (solid lines). The contour plots show 
the Si p orbital pointing in the [111] -direction between two nearest neighbors in the (211)- 
plane. Shown are the kinked partial waves (KPWs) at the three energies and the QMTO. 
The KPWs are normalized to 1, times a cubic harmonics, at the central hard sphere. The 
contours are the same in all plots. 



site and at the three energies eo, ei, and €2- In figure 2 we show this pm-QMTO, together with 
the piii-KPW at the three energies, along the hne connecting the two nearest neighbors and 
proceeding into the back-bond. 

In general, the members (labelled by R'L') of the NMTO basis set for the energy mesh 
eo,...,eAr are superpositions, 

N 

Xr'L' (r) = 2^ 't'RL (en, r) lI^Jl,r'L' ' (1) 

n=ORLeA 

of the kinked partial waves, (pjiL (s, r) , at the A^ + l points (labelled by n) of the energy mesh. In 
the present case, the L-summation is over the nine s-, p-, and d-KPWs, and the /^-summation 
is over all Si sites. Due to the localized nature of the KPWs illustrated in the figures, the 
latter summation is limited to the neighbors. The i2L-values for which we have MTOs in the 
basis set, we label active (A) , or low. Expression (|^) is the energy-quantized form of Lagrange 
interpolation, 

N N 

n=0 m.=0,yi^n 




Figure 2: Si pin-KPW for the Si spd-set plotted along the [lll]-line connecting nearest 
neighbors. The a's indicate the hard spheres at the central and the nearest-neighbor sites. 
The generating MT-potential had no repulsive potential wells at interstitial sites (E), but 
only large Si-centered wells with s = 1.7a. 



of a function of energy, (p (e) , by an A^th-degree polynomial, x^^^ i^) ■ The A^th-degree poly- 
nomial, l^'^ (e) , is substituted by a matrix with elements, L^J^^ , the function of energy, 
(j) {e) , by a Hilbert space with axes, (fiRL (e, r) , and the interpolating polynomial, x*-^"* (^) j by a 
Hilbert space with axes, X^ri\i (r) • 

As illustrated in figures 1 and 2, a kinked partial wave is basically a partial wave with a tail 
joined continuously to it with a kink at a central, so-called hard sphere of radius a^. This kink 
is seen most clearly for the lowest energy, eo- As usual, the partial wave is (pRi {e,rji)YL (fi?) , 
where the function of energy is the regular solution of the radial Schrodinger equation, 

- [r^PRi {e, r)]" = [e - vr (r) - I {I + 1) /r^] ripRi {e, r) , (2) 

for the potential-well vr (r) . The tail of the kinked partial wave is a so-called screened spherical 
wave, ipRL {£, r) , which is essentially the solution with energy e of the wave equation in the 
interstitial between the hard spheres, —Aip (e, r) = eip {s, r) , with the boundary condition that, 
independent of the energy, ipRL (e, r) go to Yl (f i?) at the central hard sphere, and to zero (with 
a kink) at all other hard spheres. It is this latter confinement, easily recognized in the plots, 
particularly at the highest energy 62, which makes the screened spherical waves, the KPWs, 
and the MTOs localized when the energy is not too high. At the same time, it makes the 
KPW have pure L-character merely at its central sphere, because outside, it is influenced by 
the hard spheres centered at the neighbors. The default value of the hard-sphere radii, qr, 
is 90 per cent of the appropriate covalent, atomic, or ionic radius. The kinked partial wave 
thus has a kink, not only at its own, but also at the neighboring hard spheres, inside which it 
essentially vanishes. 'Essentially' because the above-mentioned boundary condition only applies 
to the active components of the spherical-harmonics expansions of the screened spherical wave 
on the hard spheres. For the remaining components, in the present case the Si /- and higher 
components, as well as all components on empty (E) spheres, the screened spherical wave equals 
the corresponding partial-wave solution of Schrodinger's equation throughout the MT-sphere. 
The small bump seen in figure 1 in the lowest KPW contour along the [111] -direction is mainly 
caused by the /-character on the nearest neighbor, and so is the finite amplitude seen in figure 
2 inside the nearest hard sphere. 



3 Computational steps 

The radial Schrodinger (Dirac) equations (|2|) are integrated numerically from r = to sr. This 
yields the radial functions, ifRi (e, r) , and their phase shifts, r]Ri (e) , each of which are obtained 
by matching the logarithmic derivative of fRi (e, r) at r = sr to that of 

^Rl ^ jl i'^r) - tan rjRi (e) ni (nr) . (3) 

The radial integration must be performed for each potential well and for each I, increasing until 
all further phase shifts vanish due to dominance of the centrifugal term in (|2|) . 

The screened spherical waves are specified by a Hermitian structure matrix, whose element 
Br'l',rl (^) is essentially the radial logarithmic derivative of the L'-component in the spherical- 
harmonics expansion at the hard sphere at site R' of the screened spherical wave ipRi (Ej r) ■ 
What is known analytically, is the element 

B%L',RL (e) = Y,^7ri-'+''-'"CLL'i"Kni. {k |R - R'|) y£„ (iTHl') , 
l" 

of the bare KKR structure matrix, which specifies how the spherical wave, n; {nrR) Yl [f-R) , at 
site R is expanded around another site, R', in regular spherical waves, ji' (KrR/) Y^/ {rR') ■ Here 



Av^ = e, and the on-site terms of the bare structure matrix are defined to vanish. Screening of 
the structure matrix, 



B{e) ^] =\B^{e)^] + K-HanoBL (e) 6rr'6ll' , 



RL,R'L' 



RL,R'L' 



requires inversion of the matrix B^j^ (e) + Kcot url (s) ^rr'^lu- This can be done by fixing 
R and Hmiting R' to the 10-50 nearest sites, url (e) are the hard-sphere phase shifts for the 
active channels, 

tan aRL {e) = ji (naR) /rii {kur) , 

and, for the remaining channels, aRL (s) are the proper phase shifts, tjrl (e) . The latter channels, 
which will not have KPWs and MTOs associated with them, are said to be downfolded. With 
appropriate division into active and downfolded channels, the screened structure matrix will 
have short spatial range and no poles in the energy-range of the occupied states. 
A kinked partial wave is defined as: 

(f)RL {e, r) = [(fRi (e, tr) - (f%i (e, tr)] Yl {vr) + V'rl (e, r) , (4) 

where (e, r) is the radial solution for the central well from to s, and ip" (e, r) is the phase- 
shifted wave ^ proceeding smoothly inwards from s to the central a-sphere, where it is matched 
with a kink to the screened spherical wave ip (e, r) . The kinks of the KPW set are then given 
by the kink matrix, 

/ ^ _ Brl,r'l' (e) + K cot 7]%^ (e) Srr'6ll' 

^RL,R'L' = 7 ^ 7 ^ , (.Oj 

where 77^^ (e) is the phase shift with respect to the hard-sphere medium, 

tan (e) = tan t^rl {e) - tan url {e) . 

The rows and columns of the kink matrix run merely over active channels. In the fomalism 
above, we have for simplicity used the notation of scattering theory, which is analytical for 
e > 0, and for e < 0. Screened scattering theory with the normalization (^) is however analytical 
in a region of interest around e = 0. 

Finally, the Lagrange matrix which gives the MTO set (||) in terms of the KPW set (^), is 
given solely in terms of the values of the Green matrix, G (e) = K {e)~^, on the energy mesh 
e = eg, ei, etv- The Hamiltonian and overlap matrices in the MTO representation are given 
in terms of the same values of the Green matrix, plus the values of its first energy derivative, 
G(e). 



4 Further downfoldings 

In the valence and lowest conduction bands of Si, there are only s- and p-, but no d-electrons. 
To describe these bands, we should therefore be able to use a basis with only Si s- and p-MTOs, 
that is, with only 4 orbitals per atom. We thus let the Si s- and p-partial waves remain active, 
while the Si d- waves are now included among the passive ones, i.e. those 'folded down' into 
the tails of the screened-spherical waves in (^). The results for the bands and the pm-QMTO 
are shown in figures 3 and 4. These bands are indistinguishable from those obtained with the 
Si spd-set, on the scale of the figure, although between the energies of the mesh, the former 
bands do lie slightly above the latter. However, by making the mesh denser (increasing A^), the 
accuracy can be increased arbitrarily. The KPW of the sp-set is seen to have d-character on the 



nearest Si neighbor, and the QMTO and the KPW, particularly the one at the highest energy, 
are seen to be somewhat less localized than those for the spd-set. 

It is even possible to construct an arbitrarily accurate MTO-basis which spans merely the 
occupied orbitals, that is, which spans the valence band with a basis of merely 2 orbitals per 
atom. For tetrahedrally coordinated covalent semiconductors like Si, it is customary to take the 
valence-band orbitals as the bond-orbitals, which are the bonding linear combinations of directed 
sp'^-hybrids of orthonormal orbitals. It is, however, far simpler and more general, e.g. not limited 
to elemental semiconductors and tetrahedral structures, to take the valence-band orbitals as 
the s- and p-MTOs on every second Si atom, all partial waves on the nearest neighbors being 




Figure 3: Same as figure 1, but for the Si sp-set. 



downfolded. This corresponds to a Si^"'"Si^~ ionic picture. This QMTO-set turns out to describe 
merely the valence band, and to do so surprisingly well considering the fact that the two silicons 
are treated differently, so that the degeneracy along the XW-Iinc is, in fact, slightly broken. The 
error between the energy points is proportional to [e-j (k) — eo] [£« (k) — ei] [si (k) — £2] , exactly 
as for the basis with 4 orbitals per atom shown in figure 3, because we use QMTOs in both cases, 
but the prefactors are larger for the smaller basis: As the number of active channels decreases, 
the KPWs attain longer range and stronger energy dependence. However, by making the energy 
mesh finer, the errors of the MTO set can be made arbitrarily small. In the bottom line of figure 
5, we show the result of such a valence-band-only calculation with iV = 3 for Ge, together with 
the Pill cubic MTO (CMTO) centered on the Ge atom to the right. The accuracy of the valence 
band is superb, and the MTO is seen to spill over onto the nearest-neighbor atom(s) which were 
chosen not to have orbitals associated with them. 

Since this basis is complete for the occupied states, we may compute the density-functional 
ground-state properties in real space by taking traces, provided that we first Lowdin orthonor- 
malize the basis in real space. The sum of the one-electron energies is then computed as the 
trace of the Hamiltonian, i.e. as the sum of the energies of the orthonormal orbitals, and the 
charge density is computed as the sum of the squares of these orbitals. This is a method where 
the amount of computation increases merely linearly with the size of the system, a so-called 
order-M method. Here, M refers to the number of atoms in the system and not to the order N 
of the MTOs. This NMTO method, which generates the complete basis for the occupied states 




Figure 4: Same as figure 2, but for the Si sp-set. 



a priori, should be superior to current order-A/" methods, which either use inaccurate empirical 
tight-binding models or project onto the occupied states during the course of a large-basis-set 
calculation. 

In order to demonstrate in further detail that our method works, we consider Si in the 
diamond structure for which the valence-band Wannier functions can be taken as bond orbitals. 




Figure 5: Valence bands (solid lines) of the series CuBr - Ge, calculated with the ionic basis 
sets where the sp-MTOs are on the anion and, except for Ge, the d-MTOs are on the cation. 
The exact bands are given by dashed lines. The contour plots show the pm-MTO on the 
anion (the atom to the right). The ionicity decreases and the covalency increases from the 
top to the bottom. 



First, we orthonormalize our symmetry-breaking Si^+ Si s,px,Py,Pz QMTO set. The resulting 
Si s and Si Px orbitals are shown in the (llO)-plane in figures 6 and 7. These orthogonahzed 
orbitals are seen to remain fairly localized. Then, we transform to the four congruent sp^-hybrids 
centered nominally on every second Si atom. As figure 8 shows, such an sp^-hybrid is, in fact. 



\ 



Figure 6: Contour plot in the (llO)-plane of the orthogonahzed s-QMTO on the central 
Si atom. All partial waves on the nearest neighbors (and their lattice translations) were 
downfolded. 




Figure 7: Same as figure 6, but for the Px-QMTO. 



the bond orbital. Hence, folding all partial waves of the atoms chosen not to carry orbitals, into 
the tail of the sp^-directed orbital on one of the other Si atoms, has made that orbital look like 
a bond orbital. The reason why the figure does not show exact symmetry between the two sites 
is caused by our use of an energy mesh with only 3 points in the valence band. Making the 
energy mesh finer will generate the exact symmetry. 

Now, Sn is metal because the bonding and antibonding bands overlap. This should, however, 
not prevent our method from working for the occupied states only, because NMTOs are energy 
selective. Since the orbitals will be shaped in such a way that the basis set solves Schrodinger's 
equation exactly for the energies on the mesh, we may merely have to choose several energy 
points in the region of band overlap helow the Fermi level -and, of course, no energy points 
above. Remember that our ionic prescription does not make use of the fact that the Wannier 
functions for the valence and conductions bands are respectively bonding and antibonding. 

Since the ionic Si Si(sp) set gives the occupied states in diamond-structured Si with arbitrary 
accuracy, the same procedure with the sp-orbitals placed exclusively on the anion, and the d- 
orbitals on the cation, will of course work for any IV-IV, III-V, II- VI, and I- VII semiconductor 
and insulator. CuBr, for instance, would be thought of as an ionic compound Cu"'"Br~ with 
the closed-shell configuration Cu3d^° Br4s^4p^, and the basis should therefore have the d-MTOs 
on the Cu atoms and the s and p-MTOs on the Br atoms. This is illustrated in the upper 
line of figure 5. Being a single-site prescription, this works for CuBr in any structure. What 
we have seen is thus, that the MTO basis can be designed a priori to span the Hilbcrt space 
of the occupied states only. That is, there is one, and only one, such MTO per electron. To 
specify such a set, one would, in order to get the electron-count right, put the orbitals where 
the electrons are though to be, and leave it to the method to shape the tails of these orbitals 
in such a way that the basis solves Schrodinger's equation exactly for occupied states of the 
given static mean-field {e.g. LDA) potential. Such ionic MTO basis sets, which 'automatically' 
span the occupied -and no further- states of any hand insulator, could make density-functional 
molecular-dynamics calculations highly efficient for such systems. 




Figure 8: Same as figure 6, but for the s^?^-hybrid. For increasing iV, this orbital converges 
to the bond orbital. 



How could one imagine to treat a chemical reaction like: 2H2O 2H2 + O2 with MTO 
bases of occupied states only? For water, it is natural to use the ionic description H^O , 
according to which the orbital configuration is O p^, i.e. one would put the p orbitals on oxygen 
and fold down all partial waves centered on the hydrogens. In principle, one might stick to 
this configuration throughout the reaction, because it keeps the electron-count right. However, 
the oxygen-centered orbitals would eventually look strange and have long range, because they 
would have to separate off pieces of wave functions sitting on the hydrogens. Such oxygen 
orbitals might be more time-consuming to generate. At some stage in the reaction, it might 
therefore be appropriate to switch to configurations such as Hs | Hs |, or H Hs^ for the hydrogen 
molecule; this is analogous to our treatment of the occupied states in tetrahedrally coordinated 
Si. For O2 with the open-shell molecular configuration ppa^ ppir'^ pp-K It, we might use an 'ionic' 
configuration like: O (2; t x t| y j) O | a; t y t|) with z referring to the local z-direction of 
the molecule. 

This example immediately leads to a treatment of open-shell systems by means of spin and 
possibly orbital polarizations. We have seen that the 3rd-generation MTO method offers the 
possibility of designing single-electron bases of atom-centered localized orbitals, which span the 
wave functions in a given energy region of a given mean-field potential. These orbitals can 
even be symmetry-breaking, as in the case of diamond. Si, Ge, and Sn, without the generating 
mean field having to be so. These orbitals thus seem to have great potential in the design of 
many-electron wave functions which describe correlated electron systems in a realistic way. 

5 Conclusion 

We have solved the long-standing problem of deriving energy-independent, short-ranged or- 
bitals from scattering theory (Hubbard 1967). The present formalism contains exactly the right 
'physics and chemistry,' we feel. This should give the computational method great speed and 
accuracy, and make it a vehicle for discovery and understanding. 
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